Load libraries

library(dplyr)
library(tidyverse)
library(ggplot2)
library(ggpubr)
library(ggthemes)
library(stats) 
library(moments)
library(grid)
library(formattable)
library(gridExtra)
library(ggsignif)
library(hrbrthemes)
library(plotrix)
library(rstatix)
library(car)
library(plotly)
library(postHoc)

Load data

data <- read.csv("/Users/nuriteliash/Documents/GitHub/Sofia-varroa-project/EAG_data_all.csv")

EAG stat analysis

the tested effect (x) is the Essential oil “EO” dissolved in acetone, the “response” (y) is the electrophysiological response of the leg measured in mV. Each leg was stimulated by different stimuli in the following order:Air > Air > Air > Acetone (solvent) > 0.1 > 0.25 > 0.5 > 1 > 2.5 > 5 > Acetone (solvent) > Air

the response was normalized…. [please complete in here, by copy-pasting from the method section:)]

to test the difference in response amplitude to the different stimuli, we will use One way ANOVA. then a post-hoc Tukey test, to see which of the stimuli are significantly differnet. we followed this tutorial for ANOVA in R

Data summary

# our data contains: 
str(data)
## 'data.frame':    664 obs. of  4 variables:
##  $ leg     : Factor w/ 83 levels "VF1","VF10","VF11",..: 1 1 1 1 1 1 1 1 12 12 ...
##  $ stimuli : Factor w/ 8 levels "0.1","0.25","0.5",..: 8 1 2 3 4 5 6 7 8 1 ...
##  $ EO      : Factor w/ 9 levels "EO4309","EO4444",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ response: num  109.5 118.5 86.6 131.9 108.5 ...
table(data$stimuli, data$EO)
##              
##               EO4309 EO4444 EO4891 EO4899 EO5068 EO5663 EO5714 EO5749 EO5763
##   0.1              8     10     10     10     11      9      9      9      7
##   0.25             8     10     10     10     11      9      9      9      7
##   0.5              8     10     10     10     11      9      9      9      7
##   1                8     10     10     10     11      9      9      9      7
##   2.5              8     10     10     10     11      9      9      9      7
##   5                8     10     10     10     11      9      9      9      7
##   acet_after       8     10     10     10     11      9      9      9      7
##   acet_before      8     10     10     10     11      9      9      9      7

We have tested a total of 9 essential oils, using 83 legs. Each essential oil was tested on at least 7 different legs.

Prior to analysis, we gonna check the ANOVA assumptions:

(1) Outliers

# plot it to detect outliers by specific leg
# first sort the order of stimuli:
data$stimuli <- factor(data$stimuli,levels = c("acet_before", "0.1", "0.25", "0.5", "1", "2.5", "5","acet_after"))

box <- ggplot(data, aes(x = stimuli, y = response)) +
  geom_boxplot(aes(colour = EO)) +
  facet_wrap( ~ EO) +
  theme_linedraw() +
        ggtitle("Varroa foreleg response to different essential oils") +
        xlab("Stimuli amount (microgram)") +
        ylab("Normalized response (%)") +
        theme(axis.text.x=element_text(angle=45, hjust=1))
  
ggplotly(box, tooltip = "leg")

Remove outliered legs: “VF71”, “VF1”, “VF2”, “VF4”, “VF77”, “VF82”, “VF24”, “VF26”, “VF37”, “VF39”, “VF46”, “VF49”, “VF17”, “VF19”, “VF50”, “VF62”

data <- data %>%
  dplyr::filter(!leg %in% c("VF71", "VF1", "VF2", "VF4", "VF77", "VF82", "VF24", "VF26", "VF37", "VF39", "VF46", "VF49", "VF17", "VF19", "VF50", "VF62"))

After excluding the outliers, we proceed with the rest of the tests:

(2) Normality:

#the dependent variable should be approximately normally distributed in each cell of the design. This can be checked using the Shapiro-Wilk normality test (shapiro_test()
normality <- data %>%
  group_by(stimuli,EO) %>%
  shapiro_test(response)

normality %>%
  dplyr::filter(p<0.05)
## # A tibble: 4 x 5
##   stimuli     EO     variable statistic       p
##   <fct>       <fct>  <chr>        <dbl>   <dbl>
## 1 acet_before EO4309 response     0.761 0.0163 
## 2 0.1         EO5663 response     0.792 0.0344 
## 3 1           EO5068 response     0.783 0.0130 
## 4 acet_after  EO4891 response     0.738 0.00607

For the above four EO x stimuli combinations the response dont distribute normally (p<0.05)

(3) Homogneity of variance

# Build the linear model
model <- lapply(split(data, data$EO), function(i){
  lm(response ~ stimuli, data = i)
  })

# now you can Create a QQ plot of residuals of each EO, for example, essential oil EO4309:
ggqqplot(residuals(model$EO4309))

plot(model$EO4309, 1)

After checking the 3 assumptions, we decided to go for a non-parametric test, to test the significant difference of each stimuli concentration from the solvent.

we used a non-parametric post hoc test for multiple comparisons.

Post hoc non-parametric test for all pairwise comparisons, with Benjamini-Hochberg p-value correction

Post hoc Kruskal-Wallis rank sum test

KW <- lapply(split(data, data$EO), function(i){
  posthocKW(i$response, i$stimuli, padjust = "BH")
})

KW$EO4309$PvaluesMatrix
##             acet_before       0.1      0.25       0.5         1       2.5
## acet_before          NA        NA        NA        NA        NA        NA
## 0.1           0.9490597        NA        NA        NA        NA        NA
## 0.25          0.9490597 0.9490597        NA        NA        NA        NA
## 0.5           0.9490597 0.9490597 0.9490597        NA        NA        NA
## 1             0.9490597 0.9490597 0.9490597 0.9490597        NA        NA
## 2.5           0.9490597 0.9490597 0.9490597 0.9490597 0.9490597        NA
## 5             0.9490597 0.9490597 0.9490597 0.9490597 0.9490597 0.9490597
## acet_after    0.9490597 0.9490597 0.9490597 0.9490597 0.9490597 0.9490597
##                     5 acet_after
## acet_before        NA         NA
## 0.1                NA         NA
## 0.25               NA         NA
## 0.5                NA         NA
## 1                  NA         NA
## 2.5                NA         NA
## 5                  NA         NA
## acet_after  0.9490597         NA
KW$EO4444$PvaluesMatrix
##             acet_before       0.1      0.25       0.5         1       2.5
## acet_before          NA        NA        NA        NA        NA        NA
## 0.1           0.9132375        NA        NA        NA        NA        NA
## 0.25          0.6750838 0.6750838        NA        NA        NA        NA
## 0.5           0.6294375 0.6750838 0.8742932        NA        NA        NA
## 1             0.6289937 0.5917343 0.8742932 0.6750838        NA        NA
## 2.5           0.5917343 0.6289937 0.6750838 0.9490597 0.9490597        NA
## 5             0.2532105 0.2532105 0.6289937 0.6289937 0.6750838 0.6750838
## acet_after    0.6750838 0.6750838 0.8742932 0.8742932 0.6750838 0.9132375
##                     5 acet_after
## acet_before        NA         NA
## 0.1                NA         NA
## 0.25               NA         NA
## 0.5                NA         NA
## 1                  NA         NA
## 2.5                NA         NA
## 5                  NA         NA
## acet_after  0.6294375         NA
KW$EO4891$PvaluesMatrix
##             acet_before       0.1      0.25       0.5         1       2.5
## acet_before          NA        NA        NA        NA        NA        NA
## 0.1           0.9868483        NA        NA        NA        NA        NA
## 0.25          0.9250712 0.9868483        NA        NA        NA        NA
## 0.5           0.6312621 0.6312621 0.6324156        NA        NA        NA
## 1             0.9868483 0.9868483 0.9868483 0.6312621        NA        NA
## 2.5           0.6312621 0.6312621 0.6324156 0.9868483 0.6312621        NA
## 5             0.6312621 0.6891247 0.9250712 1.0000000 0.6312621 0.9868483
## acet_after    1.0000000 0.9868483 0.9868483 0.6312621 0.9868483 0.6312621
##                     5 acet_after
## acet_before        NA         NA
## 0.1                NA         NA
## 0.25               NA         NA
## 0.5                NA         NA
## 1                  NA         NA
## 2.5                NA         NA
## 5                  NA         NA
## acet_after  0.6312621         NA
KW$EO4899$PvaluesMatrix
##             acet_before       0.1      0.25       0.5         1       2.5
## acet_before          NA        NA        NA        NA        NA        NA
## 0.1           0.7298385        NA        NA        NA        NA        NA
## 0.25          0.6324156 0.6324156        NA        NA        NA        NA
## 0.5           0.3219956 0.3219956 0.6891247        NA        NA        NA
## 1             0.3287615 0.3219956 0.7014248 0.7927465        NA        NA
## 2.5           0.6324156 0.3219956 0.7298385 0.7298385 0.7805921        NA
## 5             0.7298385 0.7014248 0.7298385 0.6324156 0.6324156 0.6324156
## acet_after    0.6324156 0.6324156 0.7298385 0.7553550 0.7805921 0.7553550
##                     5 acet_after
## acet_before        NA         NA
## 0.1                NA         NA
## 0.25               NA         NA
## 0.5                NA         NA
## 1                  NA         NA
## 2.5                NA         NA
## 5                  NA         NA
## acet_after  0.7298385         NA
KW$EO5068$PvaluesMatrix
##             acet_before       0.1      0.25       0.5         1       2.5
## acet_before          NA        NA        NA        NA        NA        NA
## 0.1           0.8834907        NA        NA        NA        NA        NA
## 0.25          0.5393798 0.5393798        NA        NA        NA        NA
## 0.5           0.5393798 0.8946258 0.5393798        NA        NA        NA
## 1             0.5393798 0.6770675 0.5393798 0.8834907        NA        NA
## 2.5           0.5393798 0.5504130 0.5504130 0.6770675 0.8946258        NA
## 5             0.8946258 0.8834907 0.5393798 0.5504130 0.5504130 0.5393798
## acet_after    0.5393798 0.5393798 0.8946258 0.5393798 0.5917461 0.5393798
##                     5 acet_after
## acet_before        NA         NA
## 0.1                NA         NA
## 0.25               NA         NA
## 0.5                NA         NA
## 1                  NA         NA
## 2.5                NA         NA
## 5                  NA         NA
## acet_after  0.5393798         NA
KW$EO5663$PvaluesMatrix
##             acet_before       0.1      0.25       0.5         1       2.5
## acet_before          NA        NA        NA        NA        NA        NA
## 0.1           0.5917343        NA        NA        NA        NA        NA
## 0.25          0.5917343 0.8601193        NA        NA        NA        NA
## 0.5           0.5917343 0.9490597 0.9132375        NA        NA        NA
## 1             0.7936318 0.9132375 0.9490597 0.8601193        NA        NA
## 2.5           0.8601193 0.9132375 0.9132375 0.9132375 0.9132375        NA
## 5             0.9132375 0.9132375 0.8601193 0.9132375 0.9132375 0.9132375
## acet_after    0.5917343 0.9132375 0.9132375 0.9132375 0.9132375 0.8601193
##                     5 acet_after
## acet_before        NA         NA
## 0.1                NA         NA
## 0.25               NA         NA
## 0.5                NA         NA
## 1                  NA         NA
## 2.5                NA         NA
## 5                  NA         NA
## acet_after  0.8601193         NA
KW$EO5714$PvaluesMatrix
##             acet_before       0.1      0.25       0.5         1       2.5
## acet_before          NA        NA        NA        NA        NA        NA
## 0.1           0.4574500        NA        NA        NA        NA        NA
## 0.25          0.1960318 0.5256285        NA        NA        NA        NA
## 0.5           0.2223453 0.5256285 0.7332873        NA        NA        NA
## 1             0.1960318 0.3429161 0.5178863 0.5256285        NA        NA
## 2.5           0.3968159 0.6429369 0.7332873 0.8070399 0.8794139        NA
## 5             0.1960318 0.1960318 0.3429161 0.1960318 0.5178863 0.5178863
## acet_after    0.3429161 0.6429369 0.7332873 0.7194712 0.6429369 0.9490597
##                     5 acet_after
## acet_before        NA         NA
## 0.1                NA         NA
## 0.25               NA         NA
## 0.5                NA         NA
## 1                  NA         NA
## 2.5                NA         NA
## 5                  NA         NA
## acet_after  0.5178863         NA
KW$EO5749$PvaluesMatrix
##             acet_before       0.1      0.25       0.5         1       2.5
## acet_before          NA        NA        NA        NA        NA        NA
## 0.1           0.9441937        NA        NA        NA        NA        NA
## 0.25          0.9441937 0.9580911        NA        NA        NA        NA
## 0.5           0.9441937 0.9441937 0.9502984        NA        NA        NA
## 1             0.9441937 0.9502984 0.9502984 0.9502984        NA        NA
## 2.5           0.9441937 0.9441937 0.9441937 0.9502984 0.9502984        NA
## 5             0.8669914 0.8669914 0.9441937 0.9441937 0.9441937 0.9441937
## acet_after    0.8669914 0.9441937 0.9441937 0.9441937 0.9441937 0.9441937
##                     5 acet_after
## acet_before        NA         NA
## 0.1                NA         NA
## 0.25               NA         NA
## 0.5                NA         NA
## 1                  NA         NA
## 2.5                NA         NA
## 5                  NA         NA
## acet_after  0.9502984         NA
KW$EO5763$PvaluesMatrix
##             acet_before       0.1      0.25       0.5         1       2.5
## acet_before          NA        NA        NA        NA        NA        NA
## 0.1           0.9775137        NA        NA        NA        NA        NA
## 0.25          0.9775137 1.0000000        NA        NA        NA        NA
## 0.5           0.9775137 0.9775137 0.9775137        NA        NA        NA
## 1             0.9775137 0.9775137 0.9775137 1.0000000        NA        NA
## 2.5           0.9775137 0.9775137 0.9775137 0.9775137 1.0000000        NA
## 5             0.9775137 0.9775137 0.9775137 0.9775137 0.9775137 0.9775137
## acet_after    0.9775137 0.9775137 0.9775137 0.9775137 0.9775137 0.9775137
##                     5 acet_after
## acet_before        NA         NA
## 0.1                NA         NA
## 0.25               NA         NA
## 0.5                NA         NA
## 1                  NA         NA
## 2.5                NA         NA
## 5                  NA         NA
## acet_after  0.9775137         NA

Post hoc Wilcoxon test

Wilcoxon <- lapply(split(data, data$EO), function(i){
  pairwise.wilcox.test(i$response, g = i$stimuli, p.adjust.method = "BH")
})

Wilcoxon$EO4309$p.value
##            acet_before 0.1 0.25 0.5  1 2.5  5
## 0.1                  1  NA   NA  NA NA  NA NA
## 0.25                 1   1   NA  NA NA  NA NA
## 0.5                  1   1    1  NA NA  NA NA
## 1                    1   1    1   1 NA  NA NA
## 2.5                  1   1    1   1  1  NA NA
## 5                    1   1    1   1  1   1 NA
## acet_after           1   1    1   1  1   1  1
Wilcoxon$EO4444$p.value
##            acet_before       0.1      0.25       0.5        1       2.5
## 0.1          0.9708625        NA        NA        NA       NA        NA
## 0.25         0.7489510 0.7489510        NA        NA       NA        NA
## 0.5          0.7261072 0.7489510 0.9389083        NA       NA        NA
## 1            0.7261072 0.6812354 0.9389083 0.7489510       NA        NA
## 2.5          0.6812354 0.7261072 0.7489510 1.0000000 1.000000        NA
## 5            0.2447552 0.2447552 0.7261072 0.7261072 0.748951 0.7489510
## acet_after   0.7489510 0.7489510 0.9389083 0.9389083 0.748951 0.9708625
##                    5
## 0.1               NA
## 0.25              NA
## 0.5               NA
## 1                 NA
## 2.5               NA
## 5                 NA
## acet_after 0.7261072
Wilcoxon$EO4891$p.value
##            acet_before       0.1      0.25       0.5         1       2.5
## 0.1          1.0000000        NA        NA        NA        NA        NA
## 0.25         1.0000000 1.0000000        NA        NA        NA        NA
## 0.5          0.7069034 0.7069034 0.7069034        NA        NA        NA
## 1            1.0000000 1.0000000 1.0000000 0.7069034        NA        NA
## 2.5          0.7069034 0.7069034 0.7069034 1.0000000 0.7069034        NA
## 5            0.7069034 0.7645688 1.0000000 1.0000000 0.7069034 1.0000000
## acet_after   1.0000000 1.0000000 1.0000000 0.7069034 1.0000000 0.7069034
##                    5
## 0.1               NA
## 0.25              NA
## 0.5               NA
## 1                 NA
## 2.5               NA
## 5                 NA
## acet_after 0.7069034
Wilcoxon$EO4899$p.value
##            acet_before       0.1      0.25       0.5         1       2.5
## 0.1          0.7856762        NA        NA        NA        NA        NA
## 0.25         0.7069034 0.7069034        NA        NA        NA        NA
## 0.5          0.3491841 0.3491841 0.7645688        NA        NA        NA
## 1            0.3637607 0.3491841 0.7731546 0.8335142        NA        NA
## 2.5          0.7069034 0.3491841 0.7856762 0.7856762 0.8280181        NA
## 5            0.7856762 0.7731546 0.7856762 0.7069034 0.7069034 0.7069034
## acet_after   0.7069034 0.7069034 0.7856762 0.8074095 0.8280181 0.8074095
##                    5
## 0.1               NA
## 0.25              NA
## 0.5               NA
## 1                 NA
## 2.5               NA
## 5                 NA
## acet_after 0.7856762
Wilcoxon$EO5068$p.value
##            acet_before       0.1      0.25       0.5         1       2.5
## 0.1          0.9288701        NA        NA        NA        NA        NA
## 0.25         0.5946524 0.5946524        NA        NA        NA        NA
## 0.5          0.5946524 0.9314274 0.5946524        NA        NA        NA
## 1            0.5946524 0.7276018 0.5946524 0.9288701        NA        NA
## 2.5          0.5946524 0.6012341 0.6012341 0.7276018 0.9314274        NA
## 5            0.9314274 0.9288701 0.5946524 0.6012341 0.6012341 0.5946524
## acet_after   0.5946524 0.5946524 0.9314274 0.5946524 0.6429410 0.5946524
##                    5
## 0.1               NA
## 0.25              NA
## 0.5               NA
## 1                 NA
## 2.5               NA
## 5                 NA
## acet_after 0.5946524
Wilcoxon$EO5663$p.value
##            acet_before       0.1      0.25       0.5         1       2.5
## 0.1          0.6812354        NA        NA        NA        NA        NA
## 0.25         0.6812354 0.9708625        NA        NA        NA        NA
## 0.5          0.6812354 1.0000000 0.9708625        NA        NA        NA
## 1            0.9235431 0.9708625 1.0000000 0.9708625        NA        NA
## 2.5          0.9708625 0.9708625 0.9708625 0.9708625 0.9708625        NA
## 5            0.9708625 0.9708625 0.9708625 0.9708625 0.9708625 0.9708625
## acet_after   0.6812354 0.9708625 0.9708625 0.9708625 0.9708625 0.9708625
##                    5
## 0.1               NA
## 0.25              NA
## 0.5               NA
## 1                 NA
## 2.5               NA
## 5                 NA
## acet_after 0.9708625
Wilcoxon$EO5714$p.value
##            acet_before       0.1      0.25       0.5         1       2.5
## 0.1          0.5310447        NA        NA        NA        NA        NA
## 0.25         0.2121212 0.5955711        NA        NA        NA        NA
## 0.5          0.2474747 0.5955711 0.7956177        NA        NA        NA
## 1            0.2121212 0.3988604 0.5928516 0.5955711        NA        NA
## 2.5          0.4617716 0.7132867 0.7956177 0.8666846 0.9349046        NA
## 5            0.2121212 0.2121212 0.3988604 0.2121212 0.5928516 0.5928516
## acet_after   0.3988604 0.7132867 0.7956177 0.7891502 0.7132867 1.0000000
##                    5
## 0.1               NA
## 0.25              NA
## 0.5               NA
## 1                 NA
## 2.5               NA
## 5                 NA
## acet_after 0.5928516
Wilcoxon$EO5749$p.value
##            acet_before       0.1      0.25       0.5         1       2.5
## 0.1          0.9946531        NA        NA        NA        NA        NA
## 0.25         0.9946531 1.0000000        NA        NA        NA        NA
## 0.5          0.9946531 0.9946531 0.9946531        NA        NA        NA
## 1            0.9946531 0.9946531 0.9946531 0.9946531        NA        NA
## 2.5          0.9946531 0.9946531 0.9946531 0.9946531 0.9946531        NA
## 5            0.9790210 0.9790210 0.9946531 0.9946531 0.9946531 0.9946531
## acet_after   0.9790210 0.9946531 0.9946531 0.9946531 0.9946531 0.9946531
##                    5
## 0.1               NA
## 0.25              NA
## 0.5               NA
## 1                 NA
## 2.5               NA
## 5                 NA
## acet_after 0.9946531
Wilcoxon$EO5763$p.value
##            acet_before 0.1 0.25 0.5  1 2.5  5
## 0.1                  1  NA   NA  NA NA  NA NA
## 0.25                 1   1   NA  NA NA  NA NA
## 0.5                  1   1    1  NA NA  NA NA
## 1                    1   1    1   1 NA  NA NA
## 2.5                  1   1    1   1  1  NA NA
## 5                    1   1    1   1  1   1 NA
## acet_after           1   1    1   1  1   1  1

EAG plot

for some reason, i cannot add the leg ID to the hovering text in the box plot, but i can do it in the dot-plot. so in order to detect outliers, we can detect them in the first plot (the box plot), then identify their ID leg in the dot plot:)

dot <- ggplot(data, aes(x = stimuli, y = response, colour = EO, leg=leg)) +
  geom_point() + 
  facet_wrap( ~ EO) +
  theme_linedraw() +
        ggtitle("Varroa foreleg response to different essential oils") +
        xlab("Stimuli amount (microgram)") +
        ylab("Normalized response (%)") +
        theme(axis.text.x=element_text(angle=45, hjust=1)) +
        stat_summary(aes(group=EO), fun=mean, geom="line", colour="black") 

ggplotly(dot, tooltip = c("leg","response"))